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ABSTRACT 

The concentration-mass relation for dark matter-dominated halos is one of the essential results expected from 
a theory of structure formation. We present a simple prediction scheme, a cosmic emulator, for the c-M relation 
as a function of cosmological parameters for wCDM models. The emulator is constructed from 37 individual 
models, with three nested A-body gravity-only simulations carried out for each model. The mass range covered 
by the emulator is 2 • 10 12 M© < M < 10 15 M© with a corresponding redshift range of z = 0- 1. Over this range 
of mass and redshift, as well as the variation of cosmological parameters studied, the mean halo concentration 
varies from c ~ 2 to c ~ 8. The distribution of the concentration at fixed mass is Gaussian with a standard 
deviation of one-third of the mean value, almost independent of cosmology, mass, and redshift over the ranges 
probed by the simulations. We compare results from the emulator with previously derived heuristic analytic fits 
for the c-M relation, finding that they underestimate the halo concentration at high masses. Using the emulator 
to investigate the cosmology dependence of the c-M relation over the currently allowable range of values, we 
find - not surprisingly - that a% and uj m influence it considerably, but also that the dark energy equation of state 
parameter w has a substantial effect. In general, the concentration of lower-mass halos is more sensitive to changes 
in cosmological parameters as compared to cluster mass halos. 

The c-M emulator is publicly available from http : / /www . hep . anl . gov/cosmology/CosmicEmu| 

Subject headings: methods: statistical — cosmology: large-scale structure of the universe 



1. INTRODUCTION 

Over the last three decades, cosmology has made tremendous 
progress, culminating in the so-called "Standard Model of Cos- 
mology". The two main components of the Standard Model are 
a mysterious dark energy, leading to a late-time accelerated ex- 
pansion, and a dark matter component making up roughly 25% 
of the total matter-energy budget of the Universe. The evolution 
of structure in the Universe from the earliest accessible times to 
today is successfully described by a theory based on the gravi- 
tational instability - the distribution of galaxies in the Universe, 
for example, is remarkably well reproduced by this paradigm. 
Clusters and groups of galaxies are major building blocks of 
the large-scale structure and measurements of their abundance 
provide a powerful cosmological probe. Large, gravity-only A- 
body simulations have been remarkably successful in providing 
a consistent picture of the formation of the large-scale structure 
from the very early, small Gaussian density fluctuations to the 
halos, voids, and filaments we observe today. 

A surpr ising discovery from these simulations dNavarro et al.l 
fl996l[l997l) was that the dark matter-dominated halos - over a 
wide mass range typical of dwarf galaxies to massive clusters 
- share a basically universal density profile. In detail, it was 
shown that the spherically averaged density profile of relaxed 
halos formed in simulations can be described by what is now 
commonly known as the NFW (Navarro-Frenk- White) profile. 
The NFW profile is described by two parameters, the normal- 
ization and the characteristic scale radius of the halo or equiva- 
lently its (dimensionless) concentration. 

Aside from the distribution of halo masses, halo profiles are 
also of considerable interest. The profiles can be measured di- 
rectly for individual massive halos by a variety of observational 
methods, or inferred indirectly for less massive halos using sta- 
tistical lensing probes. Halo profiles are also a key input in 



halo occupation distribution (HOD) modeling of the distribu- 
tion of galaxies. NFW profiles (or minor variants thereof) are 
consistent with current observations ( Bhattachar va et ai~ll201 lb 
and, as the observations continue to improve, a corresponding 
improvement in theoretical predictions for (NFW) halo concen- 
trations as a function of cosmological parameters is needed. As 
scatter in the c-M relation is considerable, in principle, this 
would encompass knowing the actual distribution of halo con- 
centrations as a function of halo mass. 

Quantitative predictions for the c-M relation from a first 
principles analytic approach are difficult to obtain, due to the 
highly nonlinear dynamics involved in the formation of halos. 
Accurate predictions can only be obtained from computation- 
ally expensive, high-resolution simulations. These simulations 
need to cover large volumes in order to yield good statistics, 
especially in the cluster mass regime, as well as high force res- 
olution to reliably resolve the halo profiles. In recent years, 
the focus has therefore been on generating predictions for one 
cosmology around the best-fit WMAP (Wi lkinson Microwave 
Anisotropy Probe) results of that time (e.g..|Duffv et ai~l l2008. 
Bhat tacharva et al. 11201 UlPrada et al. 1120120 . The fitting func- 
tions so generated cannot be extended beyond the cosmolog- 
ical model they have been tuned for. Heuristic models tha t 
aim to extend this re ach, e.g. those by iBullock et al. I (12 001) 
and Eke et al.l d2001l) and improvements thereof ( Maccio et al. I 
20081 ) do not lead to the desired accuracy, as discussed in 
Duffv^etaljj2008). Other discussions of this issue can be 



found in iGaoetal.l d2008l) . lHavashi & White I d2008l) . and 
IZhao et al. I d2009l) . 

In order to overcome the many shortcomings of fitting func- 
tions as a general approach in cosmology, we have recently de- 
veloped the "Cosmic Calibration Framework" (CCF) to pr ovide 

accurate prediction schem es for cosmological observables ( He itmann et a, 
2006; lHabib et al.l 120071) . The aim of the CCF is to build 
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codes that act as very fast - basically instantaneous - prediction 
tools for large scal e structure observables such as the nonlinear 
powe r spectrum dHeitmann et al.1 1201 (j, 120091; lLawrence et al.l 
l201Ch . mass functions for different halo definitions, or the 
concentration-mass relation - as discussed here. Predicting 
these observables requires running a number of high-performance 
simulations to reliably resolve the nonlinear regime of structure 
formation. The CCF provides a powerful way to build preci- 
sion prediction tools from a limited number of computationally 
expensive simulations. 

At the heart of the CCF lies a sophisticated sampling scheme 
that provides an optimal sampling strategy for the cosmolog- 
ical models to be simulated (we use orthogonal array-based 
Latin hypercube as well as symmetric Latin hypercube designs; 
an introduction to the general sampling strategy is provided 
in ISantner et al.l 120031) . an optimal representation to translate 
the measurements from the simulations into functions that can 
be easily interpolated (a principal component basis turns out to 
be an efficient representation), and finally a very accurate inter- 
polation scheme (our choice here is Gaussian process model- 
ing). 

The CCF was first introduced in Heit mann et al.l (2006) an d a 
more d etailed description and examples are provided in Habib et ail 
(12001) . In a series of three papers (Coyote Universe I-III) we 
developed an emulator for the matter power spectrum for a five 
dimensional parameter space covering 9 = {ujb,u) m ,n s ,w,cr%}. 
This emulator provides predictions for the power spectrum for 
wCDM cosmologies out to k ~ 1 Mpc -1 at th e 1% accuracy 
level f or a redshift range of < z < 1 . In ISchneider et al.l 
(2008) the work was extended to derive an approximate sta- 
tistical model for the sample variance distr ibution of the non- 
linear matter power spectrum. lEifleri d2012l) used the emulator 
to generate a weak lensing prediction code to calculate various 
second-order cosmic shear statistics, e.g., shear power spec- 
trum, shear-shear correlation function, ring statistics and Com- 
plete Orthogonal Set of EB-mode Integrals (COSEBIs). 

The focus of this paper is the development of an emulator 
for the concentration-mass relation for wCD M cosmologies. 
We us e the same base set of simulations as in lLawrence et al.l 
(2010), consisting of 37 cosmological models and a single 
1300 Mpc volume, high-resolution simulation for each model. 
This simulation set is augmented here with a set of new, higher 
resolution simulations. These simulations cover smaller vol- 
umes (a 360 Mpc and a 180 Mpc simulation for each model) 
to obtain good statistics over a large range of halo masses. For 
each model we measure the best- fit concentration-mass (c-M) 
relation, assuming a simple power law form. The fits lay the 
foundation for building the emulator that provides predictions 
for the c-M relation within the wCDM parameter space cov- 
ered by the original simulations. In redshift, the emulator cov- 
ers the range between z = and z = 1 . We provide a fast code 
that delivers the mean c-M relations for wCDM cosmologies 
to good accurac^Q. 

As is well-known, the c-M relation has considerable scatter 
and, in principle, it is not obvious that this scatter should have a 
simple form, and what it s cosmological dependence might be. 
However, as discussed in Bhattacharya et al. (201 1), the scat- 
ter has a simple Gaussian form in wCDM models, and more- 
over, even though the mean c-M relation is clearly cosmology- 
dependent, as is the associated concentration variance, ct^(M), 
the ratio of <r c (M) to the mean concentration is close to 1/3, 

* |http : / / www. hep . anl . gov/ cosmology/ CosmicEmu| 



independent of cosmology, mass, or redshift. This means that 
once an emulator for the c-M relation is in hand, the concen- 
tration standard deviation is given automatically by a simple 
relation. 

The paper is organized as follows. After a brief outline of 
the halo concentration measurements from the simulations, we 
describe the cosmological model space and the simulation suite 
used to build the emulator. In Section [3] we also discuss the 
generation of the smooth prediction for the concentration-mass 
relation for each model that underlies the interpolation scheme 
for building the emulator. We give a brief description on how to 
build the emulator in Section|4]and show some examples from 
the working emulator and test results verifying its accuracy. We 
also compare our results to currently used fitting formulae and 
investigate the cosmology dependence of the c-M relation in 
some detail. Finally, we provide a conclusion and outlook in 
Section|5] 

2. CONCENTRATION-MASS RELATION 

We study the concentration-mass relation in the regime of 
bright galaxies to clusters of galaxies, spanning halo mass 
ranges between 2 ■ 10 12 M© to 10 15 M Q , while varying wCDM 
cosmological parameters. A detailed description on how to 
measure halo concentrations from simulations and a discussio n 
of possible systematics is given in Bhattachary a et al."1 d201 lb . 
We follow the same approach in this paper and give here a brief 
summary of the main steps in measuring the c-M relation in 
our simulations. 

As a first step, we identi fy halos using a fast p arallel friends- 
of-friends (FOF) finder dWoodring et al. I [201 ll) with linking 
length b = 0.2. Once a halo is found, we define its center via a 
density maximum criteria - the location of the particle with the 
maximum number of neighbors. This definition of the halo cen- 
ter is very close to that given by the halo's potential minimum. 
Given a halo center, we grow spheres around it and compute 
the mass in radial bins. Note that even though an FOF finder 
is used, the actual halo mass is defined by a spherical overden- 
sity method, consistent with what is done in observations. (For 
discu ssion s on halo mass, se e, e.g., I White 1 1200 U iLukic et al. I 
20091 and lMoreetaT~ll2011h . The NFW form for the spheri- 
cally averaged halo profile is a function of two parameters, one 
of which is constrained by the halo mass. Here we fit the mass 
profile using both total halo mass and concentration as free vari- 
ables. Although the mass could be measured independently of 
the concentration, the joint analysis is potentially less sensitive 
to fitting bias. 

We write the NFW profile as 

1 \ $ Petit 

(r/r s ){\ + r/r s y 

where 8 is a characteristic dimensionless density, and r s is the 
scale radius of the NFW profile. The concentration of a halo is 
defined as ca = r& / r s , where A is the overdensity with respect 
to the critical density of the Universe, p cl -i t = 3H 2 /8irG, and 
r& is the radius at which the enclosed mass, Ma, equals the 
volume of the sphere times the density Ap cl -i t . We compute 
concentrations corresponding to A = 200, corresponding in turn 

tOC200 =^200/r s - 

The mass enclosed within a radius r for an NFW halo is given 

by 

M(< r) = m(c A r/R2oo)/m(c2oo)M 2 oo, (2) 
where m(y) = ln(l +y) — y/(l +y). The mass in a radial bin is 
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Table 1 



# 



The parameters 



FOR THE 



37 + 1 MODELS WHICH 
—W <J% 



DEFINE 

h 



THE SAMPLE 
# 



SPACE 



See text 



FOR FURTHER DETAILS. 

n s —w (78 



M000 
M001 
M002 
M003 
M004 
M005 
M006 
M007 
M008 
M009 
M010 
M011 
M012 
MOO 
M014 
M015 
M016 
M017 
M018 



0.1296 
0.1539 
0.1460 
0.1324 
0.1381 
0.1358 
0.1516 
0.1268 
0.1448 
0.1392 
0.1403 
0.1437 
0.1223 
0.1482 
0.1471 
0.1415 
0.1245 
0.1426 
0.1313 



0.0224 
0.0231 
0.0227 
0.0235 
0.0227 
0.0216 
0.0229 
0.0223 
0.0223 
0.0234 
0.0218 
0.0234 
0.0225 
0.0221 
0.0233 
0.0230 
0.0218 
0.0215 
0.0216 



0.9700 
0.9468 
0.8952 
0.9984 
0.9339 
0.9726 
0.9145 
0.9210 
0.9855 
0.9790 
0.8565 
0.8823 
1.0048 
0.9597 
1.0306 
1.0177 
0.9403 
0.9274 
0.8887 



1.000 
0.816 
0.758 
0.874 
1.087 
1.242 
1.223 
0.700 
1.203 
0.739 
0.990 
1.126 
0.971 
0.855 
1.010 
1.281 
1.145 
0.893 
1.029 



0.8000 
0.8161 
0.8548 
0.8484 
0.7000 
0.8226 
0.6705 
0.7474 
0.8090 
0.6692 
0.7556 
0.7276 
0.6271 
0.6508 
0.7075 
0.7692 
0.7437 
0.6865 
0.6440 



0.7200 
0.5977 
0.5970 
0.6763 
0.7204 
0.7669 
0.7040 
0.6189 
0.7218 
0.6127 
0.6695 
0.7177 
0.7396 
0.6107 
0.6688 
0.7737 
0.7929 
0.6305 
0.7136 



M019 0. 

M020 0. 

M021 0. 

M022 0. 

M023 0. 

M024 0. 

M025 0. 

M026 0. 

M027 0. 

M028 0. 

M029 0. 

M030 0. 

M031 0. 

M032 0. 

M033 0. 

M034 0. 

M035 0. 

M036 0. 

M037 0. 



1279 
1290 
1335 
1505 
1211 
1302 
1494 
1347 
1369 
1527 
1256 
1234 
1550 
1200 
1399 
1497 
1485 
1216 
1495 



0.0232 
0.0220 
0.0221 
0.0225 
0.0220 
0.0226 
0.0217 
0.0232 
0.0224 
0.0222 
0.0228 
0.0230 
0.0219 
0.0229 
0.0225 
0.0227 
0.0221 
0.0233 
0.0228 



0.8629 
1.0242 
1.0371 
1.0500 
0.9016 
0.9532 
1.0113 
0.9081 
0.8500 
0.8694 
1.0435 
0.8758 
0.9919 
0.9661 
1.0407 
0.9239 
0.9604 
0.9387 
1.0233 



1.184 
0.797 
1.165 
1.107 
1.261 
1.300 
0.719 
0.952 
0.836 
0.932 
0.913 
0.777 
1.068 
1.048 
1.147 
1.000 
0.853 
0.706 
1.294 



0.6159 
0.7972 
0.6563 
0.7678 
0.6664 
0.6644 
0.7398 
0.7995 
0.7111 
0.8068 
0.7087 
0.6739 
0.7041 
0.7556 
0.8645 
0.8734 
0.8822 
0.8911 
0.9000 



0.8120 
0.6442 
0.7601 
0.6736 
0.8694 
0.8380 
0.5724 
0.6931 
0.6387 
0.6189 
0.7067 
0.6626 
0.6394 
0.7901 
0.7286 
0.6510 
0.6100 
0.6421 
0.7313 



then 



M i = M{<r,)-M{<r i - 1 ). 



(3) 



We then fit Eq.[3]to the mass contained in the radial bins of each 
halo, by minimizing the associated value of \ 2 as 



X 



E 



(Mf m -Mj) 
(M,"'"') 2 /«/ 



(4) 



where the sum is over the radial bins, n, is the number of parti- 
cles in a radial bin, Mf m is the mass in bin i calculated from the 
simulations and M, is the mass calculated assuming the NFW 
profile. The advantage of fitting mass in radial bins rather than 
the density is that the bin center does not have to be specified. 
Note that we explicitly account for the finite number of parti- 
cles in a bin. This leads to a slightly larger error in the profile 
fitting but minimizes any possible bias due to the finite number 
of particles, especially near the halo center. 

We fit for two parameters - the normalization of the profile 
and the concentration. Halo profiles are fitted in the radial range 
of approximately (0.1 - l)R V j r - This choice is motivated partly 
by the observations of concentrations that typicall y exclude the 
centra l region of clusters (e.g., observations by lOguri et al. I 
l20Tll) . More significantly, however, this excludes the central 
core which is sensitive to the effects of baryonic physics and 
numerical errors arising from limitations in both mass and force 
resolution. iDuffv et al. I d2010l) have shown that, at r < 0.1/? V! >, 
cluster halo profiles are potentially sensitive to the impact of 
baryons with the profiles being affected at r = 0.05R V j r by as 
much as a factor of 2. 

The c-M relation is calculated by weighing the individual 
concentrations by the halo mass, 



C(M) : 

where the sum is over the number, Nj, 
bin. The mass of the bin is given by 



(5) 



of the halos in a mass 



(6) 



The error on c(M) is the mass-weighted error on the individual 
fits plus the Poisson error due to the finite number of halos in 
an individual bin added in quadrature, 



Ac(M) : 



(7) 



V EiMi J Ni 

where Ac,- is the individual concentration error for each halo. 
The first term dominates towards the lower mass end where the 
individual halos have smaller number of particles and the sec- 
ond term dominates towards the higher mass end, where there 
are fewer halos to average over. 

3. COSMOLOGICAL MODELS AND SIMULATION SETS 

We now describe the cosmological model space covered by 
our prediction scheme and the simulations used to construct it. 
The emulator is based on 37 cosmological models spanning the 
class of wCDM cosmologies. We allow for variations of the 
following five parameters: 

0= { w,,, w,„, n. 5 ,w, cr 8 }• (8) 

The 37 models are chosen to lie within the ranges: 

0.0215 <LJ b < 0.0235, 
0.120 <ui m < 0.155, 

0.85 <n s < 1.05, (9) 
-1.30 < w<-0.70, 
0.616 < cr 8 <0.9, 

which are picked based on curre nt constraints from CMB mea- 

surem ents dKomatsu et al.l20l"Tl) . Following the approach in lLawrence et ; 
(2010) we lock the value of the Hubble parameter h to the best- 
fit value for each model, given the measurement of the distance 
to the surface of last scattering. The values for h then range 
from 0.55 < h < 0.85. In addition to the 37 models, we run 
one ACDM model (M000 in Table [TJ which is not used to build 
the emulator. Instead we use this model as a control for testing 
the accuracy of the emulator. All 37+1 models are specified in 
detail in Table [TJ 
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FIG. 1. — Concentration-mass relations for 37 wCDM cosmologies. The blue points show the measurements from the three simulations per model while the red 
lines show the best-fit power law for each measurement. In each subplot we show the results for z = and 2=1 (upper and lower curve respectively). We find the 
best-fit power law separately for both redshifts. Models with low values for erg in general also exhibit lower c—M relations (e.g. M012, M018, M019). The fits 
shown here are the foundation for building the emulator described in Section|4] 



The specific model se lection process is described at length 
in lHeitmann et al.l d2009h . In summar y, it is based on Symmet- 
ric Latin Hypercube (SLH) sampling dLi & Yej20 00): this sam- 
pling strategy provides a scheme that guarantees good coverage 
of the parameter hypercube. In our specific case we choose an 
SLH design that has good space filling properties in the case 
of two-dimensional projections in parameter space. In other 
words, if any two parameters are displayed in a plane, the plane 
will b e well covered by simulation points. Heitmann et al.l 
( 120091) provide an extensive discussion regarding optimal de- 
sign choices and we refer the interested reader to that paper. 

The emulator developed here is valid between < z < 1 



Table 2 

box sizes, particle numbers, and mass resolution. 



Length [Mpc] 




Force res. [kpc] 


m p 


[M ] 


1300 


1024 3 


50 


5.7- 


\Q U LU m 


365 


512 3 


10 


1.0- 




180 


512 3 


10 


1.2- 


10 10 w m 



and covers a halo mass range from 2 • 10 12 M Q to 10 15 Mq. 
We use different box sizes to cover different mass ranges with 
sufficient statistics. A summary of the different simulation 
sizes is given in Table [2] All simul ations were ca rried out 
with the Tre ePM code GADGET- 2 jSpringe]||2005al) . In pre- 
vious work dBhattacharya et ai~l 1201 ll) . we have shown that 
results from GADGET- 2 simulation s and those with HACC 
dHabib et al. 2009: [Pope et al. 11201 Ol) produce completely con- 
sistent results. Results from a recent cluster re-simulation cam- 
paig n dWu et al. 20121) are a lso in good agreement with those 
of Bh attacharya et al. I (1201 ll) . 

One set of simulatio ns is from the origina l Coyote Universe 
suite as described in lLawrence et al. | d2010l) . This set of runs 
evolves 1024 3 particles in (1300 Mpc) 3 volumes. In addition, 
we run one realization each per model with 512 3 particles with 
a 10 kpc force resolution in a 365 Mpc box and a 180 Mpc 
box. A summary of the simulation sets including force and 
mass resolution is given in Table 

We combine the simulation results from the three boxes for 
each model to obtain measurements spanning the desired mass 
range. While some models (in particular those with high values 
of erg) have clusters at even higher masses, the statistics beyond 
10 15 M Q are insufficient and we exclude those measurements. 
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Fig. 2. — Ratio of the standard deviation of the concentration to the mean 
concentration as a function of mass for all 37 cosmologies at z = 0. While both 
the standard deviation and the mean concentration are functions of cosmology 
and redshift, their ratio is essentially invariant, and is approximately 1/3. The 
distribution of the concentration around the mean is well-fit by a Gaussian 
distribution I Bhattacharya et al. 201 1). 

This is also done in order to avoid extrapolations for models 
where no data points at high masses exist. In order to build an 
emulator, for each model we have to provide a prediction for 
the c-M relation for the same mass range. This ensures that 
we can provide a consistent set of measurements for the final 
interpolation process between different models. From the sim- 
ulation results, we determine for each of the 37+1 cosmologies 
the best-fit c-M relation by simply finding the best-fit power 
law for each model at two redshifts, z = and z = 1 . The re- 
sults for the 37 models underlying the emulator are shown in 
Fig. Q] The blue points show the simulation results while the 
red curves show the best-fit power law for each model. The 
upper curves in each plot are obtained at redshift z = and the 
lower curves at z = 1 ■ The concentration values range between 
c ~ 2 and c ~ 8. As expected, we find that models with low 
values of cr 8 (e.g., M012, M018, M019 with cj 8 < 0.65) have 
depressed c-M relations. We will return to the cosmology de- 
pendence of the c-M relation in Section l4~2l after constructing 
the emulator, which will allow us to carry out a comprehensive 
sensitivity analysis. We reiterate that the fits shown in Fig.[T]are 
the basis for building the emulator; this procedure is discussed 
in the next section. 



Finally we turn to a discussion of the intrinsic scatter in the 
c-M relation. As mentioned earlier, the distribution of con- 
centrations at any given halo mass is Gaussian, and the ratio 
of a c {M) to the mean concentration is an ap proximate invariant 
for w CDM models, with a value of ~ 1/3 (Bhattach arya et al. I 
1201 ll) . independent of redshift and halo mass. This behavior is 
exhibited in Fig. [2] where the ratio is computed for all 37 cos- 
mologies as a function of halo mass, at z = 0. Thus, given the 
c-M relation from the emulator, the standard deviation at each 
mass bin can be trivially estimated by multiplying the returned 
concentration value by 1/3. 

4. EMULATOR FOR THE CONCENTRATION-MASS RELATION 
4. 1 . Building the emulator 
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FIG. 3. — First four PC basis vectors, </>,-. Absolute values are used to show 
the dynamic range on a logarithmic scale. Only the first three basis vectors are 
actually used in the emulation; any others contribute on a scale many orders of 
magnitude smaller. 



In this section, we briefly outline the process for build- 
in g the c-M emulator. We follow the procedure explained 
in Heitmann et al. (2009) and refer the reader to this pape r for 
more complete details. The focus o flHeitmann et alj d2009) was 
on modeling the matter power spectrum rather than the c—M 
relation, however, the process is essentially unchanged. Start- 
ing with the design of 37 models given in Table Q] we measure 
the c-M relation for each cosmology at z = and z = 1 and 
fit these with a power law as described in Section [3] to obtain 
a smooth functional form. First, for every mass bin, the global 
mean value is subtracted, and then via a simple rescaling, the 
concentrations are normalized to have unit variance. This pro- 
duces a zero-mean, unit variance dataset spanning the 37 cos- 
mologies. To reduce the dimensionality of the problem, these 
normalized functions are then decomposed into principal com- 
ponent (PC) basis functions and only the most significant com- 
ponents are kept. The idea is to apply the interpolation method 
of choice (Gaussian process modeling in our case) to the coef- 
ficie nts of the basis func tions, rather than to the raw data itself 
(see lHeitmann et ai]|2009l for details). Figure [3] shows that we 
only need three PCs to successfully capture the behavior of the 
c—M relation since the shape of the relationship remains fairly 
simple across this set of cosmologies. As explained above, the 
emulator actually returns the weight on each PC basis function 
and these can be combined together to give the new c-M rela- 
tion. We model the error in the projection to the PC basis with 
an additional hyperparameter X p that can be tuned to represent 
the level of noise in the data. 



A Gaussian process is then used to interpolate between the 
model results; this means that the c-M relation for a new cos- 
mology is actually a function drawn from a unit normal distri- 
bution. The covariance matrix describes the 'distance' between 
the new model and the set of known models as given by the co- 
variance function. The full covariance matrix, S, is composed 
of one £/ for each PC, arranged along the diagonal elements 
such that: £ = diag(E\...T, n ) for n PCs. Each element of S; is 
given by: 

^■,,j = ^t[pT k ~ ejt) \ do) 

k=l 
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where Q\ represents the cosmological parameters and the i and 
j indices run over the number of models spanning the design 
space (in this case i, j = 1 - 37), the I index runs over the num- 
ber of PCs and the k index runs over the number of cosmo- 
logical parameters. The hyperparameters, Xi,pu,X P , are set by 
exploring the likelihood surface, which is done with a Markov 
chain Monte Carlo analysis, but any other algorithm that lo- 
cates the maximum likelihood of a multidimensional surface 
could also be used. The complet e expression for the po sterior 
can be found in Equation B17 of lHeitm ann et aQ d2009l) . This 
conditions the Gaussian process to the design of the 37 mod- 
els and ensures that the hyperparameters correctly capture the 
complexity of the surface, because they control the fit of the 
interpolating functions to the data. 

After conditioning the GP for the best-fitting hyperparame- 
ters, the emulator is ready to predict the c-M relation for a 
different cosmology. The prediction involves re-calculating the 
covariance matrix between the new parameters and the design 
and this locates the new parameters within the design space. 
This process is quite fast, and can be repeated each time a new 
cosmology is needed with little computational cost. 

The results at intermediate redshifts (0 < z < 1) are produced 
with a simple linear interpolation. This remains fairly accurate 
because the change in the c-M relation with redshift is largely 
a simple shift in amplitude. 

4.2. Testing the Emulator 

The accuracy of the emulator is determined using two meth- 
ods: 1) we compare the performance of the emulator against a 
model not included in the original design and 2) we remove one 
of the models from the design and rebuild the emulator based 
on the remaining 36 models in what is known as a holdout test. 
In this section, we perform both of these tests to demonstrate 
the accuracy of the c-M emulator. 

We withheld one model (M000) with a ACDM concordance 
cosmology from the set of 37 models when building the c-M 
emulator. Figure|4]shows the comparison between the emulator 
prediction for this cosmology against the direct simulation re- 
sults from three different box sizes at z = and z = 1 . The hashed 
region covers the 1-sigma boundary around the mean. The em- 
ulator predictions are consistent with the A^-body c-M relations 
well within the errors on the measurements. In comparison with 
the smoothed fit for model M000, derived from the same power 
law fitting procedure used on the set of 37 cosmologies, we find 
that at z = 0, the emulator is essentially perfect at the high mass 
end and accurate to at least 3.25% for low mass halos. For 
z = 1 the error is somewhat worse, mainly due to the limited 
halo statistics for building the emulator, especially for the low- 
ers models. At low masses the predictions are accurate at the 
2% level and degrade to 9% inaccuracy at the highest masses 
considered. All of these values are well within what may be 
considered to be the nominal uncertainty in determining con- 
centrations from simulations dBhattacharva et aL~ll201 ll) For 
most of the range of halo mass considered, the accuracy of the 
emulator outperforms any other prediction scheme available, 
especially considering the large model space covered here. 

In Fig. [3] we show estimates of the emulator error by per- 
forming a holdout test. In such a test one model is kept aside 
and a new emulator is built, based on the remaining 36 models. 
The new emulator is used to predict the c-M relation for the 
held-out model. Since the numerical result ('truth') is known 
for that model, we can measure the emulator prediction error. 




13.5 
log10(MJ 
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FIG. 4. — Predictions from the c-M emulator at z = (blue, solid) and 
Z = 1 (red, solid) for measurements from iV-body simulations for the M000 
cosmology. The dashed lines show the best-fit power law describing the N- 
body results and the hashed region shows the expected variation in the c—M 
relation from the mean (solid line). Note that this set of simulations was not 
used to build the emulator. 



One shortcoming of this method - in particular if only a very 
small number of simulations is available as is the case here - is 
that by removing one model, the quality of the emulator is de- 
graded. Therefore, the error estimate for the emulator obtained 
this way can be considered to be a conservative upper bound. 

We have chosen to exclude only models M004, M008, MO 13, 
MO 16, M020 and M026, because these are located relatively 
close to the center of the design. Removing a model that de- 
fines one of the edges of the design would greatly reduce the 
performance of the emulator, since the GP would be extrapolat- 
ing for a missing model that is now outside of the design range. 
The comparison is made with respect to the smoothed A^-body 
result that was used to construct the full emulator, not the raw 
concentration measurements from the simulation. At most, the 
emulator deviates by 3.3% from the simulation results at z = 
and this rises to 15% at z = 1. This is because the error on the 
raw measurements increases with redshift as the sample size of 
halos decreases, particularly for low-cg models. 

4.3. Comparison with other c—M Predictions 

We now compare the results obtain ed from the emula t or wit h 
those from the models presented inlBhattacharva et al. I d201 j]). 
Bullock et al71d200ll).lDuffv et al. 1 (2008) and lPrada et al. I(|2012|) . 



The iBullock et al. I (2001) model was intended to correct the 
redshift dependence of the original NFW model, which was 
claimed to overpredict the concentration of high redshift (z > 1 ) 
halos. We perform our comparisions against the most re- 
cent version of the model t hat incorporates co rrections from 
IMaccio et al.1 d2008lR The IBullock et all (120011) model con- 
tains two free parameters K = 3.85 and F = 0-01- Ne wer values 
of K and F were obtained in IMaccio et al. I (12008) by fitting 
this model to A^-body simulations using cosmological parame- 
ters corresponding to the first, third and fifth WMAP data re- 
leases. Figure [6] shows the ratio of the IBullock et al. I (12001 1) 
model to our emulator at z = fo r two cosmologies, M00 and 
WMAP7 dKomatsuet al.l 1201 ll) ; note that the IBullock et al. I 
d2001l) model has only been tested with ACDM and SCDM 
cosmologies. These two models are certainly consistent at low 

2 available from physics . uci . eclu/~bullock/CVIR 
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FIG. 5. — Holdout tests for the c-M emulator at z = (top) and z = 1 
(bottom). In both plots, the labelled model has been removed from the design 
and an emulator is rebuilt on the reduced design. We then take the ratio of the 
smoothed W-body result and the prediction from the new emulator to check the 
accuracy of the full emulator made with the original design. 



halo masses, within the expected error of the emulator, but a 
substantial discrepancy occurs at cluster-size d halos, even with 
the updated version of iMaccio et al. I d2008l) . This occurs be- 
cause the model contains free parameters that need to be tuned 
to a particular cosmology with A^-body simulations. However, 
the iBullock et al. I (120011) model is able to reach much lower 
halo masses, M < 10 10 M Q , than our emulator because it is 
calibrated to higher mass resolution A^-body simulations. 

At z = 1, the public code used for the Bullock/Maccio model 
fails to compute the concentration across the full range of halo 
masses because of difficulties at low erg. We therefore show 
only results for a limited mass range. The discrepancy here is 
much larger than for z = 0, with a concentration underestimation 
of greater than 20%. 

More recently. iDuffv et al. I d2008l) proposed a new c-M re- 
lation with a power-law relationship between the halo mass 
and the concentration, as extracted from a series of high res- 
olution, small to medium volume A^-body simulations with a 
WMAP5 cosmology. Results from the c-M emulator are con- 
sistent with t heir predictions to within ~ 10% at z = 0, as are 
the results in Bhattacha rya et al."1 d201 ll) . There is a slight de- 
viation at cluster sized halos; our c-M emulator is based on 
larger volume simulations, and is therefore able to provide a 
more complete sample of massive halos and reduced shot noise 



FIG. 6. — Comparison of the emulator - taken as the reference - with 
other models for the c — M relation at z = (uppe r panel) and z = 1 (lower 
panel). The lBullock et al. K 2001),Macci6 et afl j |2008h c-M relation is shown 
in blue for both M000 ( solid) and WMAP7 (dashed) cosmologies. The 
Bhattachar ya et al. 1 1201 lh fit (note that this fit wa s derived only for M000, 
not for general cosmologi cal models) i s in green, me lPradaetaT1l20H) c-M 
relation in cyan, and the Duffy et aQ 420081) c — M relation in red for M000 
(solid) and WMAP5 (dashed) cosmologies. In addition, the pink line shows 
the ratio of the power-law fit for M000 to the emulator prediction (As shown 
in Fig. f4]the agreement is very good, with less than 3% deviation over most 
of the mass range). The l ower panel shows the results for z - 1 ■ The publicly 
available code for the Bull ock et al. I <2001hlMacci6 et al. 1 0008) fit does not 
work seamlessly over the full mass range so we do not show results for the very 
high mass end here. See the text for further discussion of this set of results. 



at th e high mass end. T he agreement between the emulator and 
the IDuffv et al. I (120081) c-M relation improves to ~ 5-6% at 
z = 1 for the WMAP5 cosmology. 



In Figure [6] we also show the ratio between our emulator 
an d the c-M relation as determined by the model discussed 
in iPrada et aL~l (120121) . which is itself based on a number of 
A^-body simulations. There is a ~ 20% discrepancy at z = 
(~ 40% at z = 1) for lower halo masses. This increases dra- 
mati cally for clu s ter siz ed haloes at both redshifts, since un- 
like IPrada et al. I (120121) . we do not observe an upturn in the 
c-M relation, where their concentration increases with halo 
mass. One should note that the methods for measuring the halo 
concentration are different in our two cas es - we use a finite- 
range profile-fittin g method as discu ssed in [Bhattach arya et al. I 
d201ll) . whereas Prada et aT] (1201 2|) use a two-point ratio method. 
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FIG. 7. — Sensitivity of the c—M relation to the five cosmological parameters varied in the emulator design, f2 m /i 2 , Sli,/i 2 . — w, n s and crs, at z = 0. We vary each 
cosmological parameter individually for each panel and have binned the range into five intervals, which are coloured from light to dark as the value of the parameter 
increases. From each of these, we subtract the c—M relation for the model corresponding to the midpoint of the design space, c(M)o- To guide the eye, we also plot 
the median c(M) - c(M)o of each bin in parameter space. 



Exectations for discrepancies between profile fitting and their 
pa rticular ratio method are f urther discussed in the appendix 
of lBhattacharva et al. I (|201 ll) . 

We also note that there is a good agreement to within ~5% 
betw een our emulator and the c—M relat ions measured bylNeto et al. 
(2007) from the Millennium simulation Springel (20053). The 
cosmology of the Millennium simulation does not quite fall 
within the range of our emulator (oV"* 2 = 0.024 and h = 0.73), 
and to facilitate this comparsion, we have adoped values as 
close to these as possible that still lie within our parameter 
space (u b h 2 = 0.0235 and h = 0.719). 

Lastly, we compar e the emulator prediction with the fitting 
function derived in Bhattac harva et al. I d2011l) for the M000 
cosmology. Before doing so, we provide some necessary back- 
ground. First, the r edshif t dependence in the fitting form in 
Bhattach arva et al."1 d201ll) is handled differ ently than in the 
current paper. In Bhattacha rva et al71 d201 lb . the aim was to 
find a global power-law fit that encompasses all redshifts con- 
sidered (between z = and z = 2) at once. Therefore, the fit for 
each redshift is not expected to be perfect. In the current paper 
we follow a different path: since we do not provide a single for- 
mula for the c-M relation but rather a simple numerical code, 
we can generate the best-fit power-law model for each redshift 
separately and then simply interpolate between the redshifts. 
This produces a more accurate answer at each redshift at the 
minimal cost of running a fast code for every c-M prediction 
instead of using one fitting formula. 

Second, for the high-mass range, Bhattachar va et al."1 d201 ll) 
used high er force resolution simulati ons. As shown in the Ap- 
pendix of Bhattac harva et al.~l d201 ll) . the concentrations from 
the Coyote runs are slightly lower at high masses (at the 5% 
level) compared to higher-resolution simulations. Since it is 
not clear if this effect is independent of cosmology (most likely 
for lower erg simulations the effect will be smaller) we decided 
to not attempt to correct the concentration measures in this pa- 
per for the Coyote runs. Therefore, the uncertainty for the 
high mass concentrations from the emulator predictions will 
be slightly higher and one expects the predictions to be biased 
slightly low. Considering the overall scatter and uncertainty in 
the c-M relation, this small effect is unlikely to be significant. 
One should note, however, that due to this suppression, the ra- 
tio of cr c (M) to the mean concentration, which we quote at the 
nominal value of 1/3, could be slightly smaller at higher masses. 

Keeping these c aveats in mind, we now t urn to the compar- 
ison of the fit by Bhatt acharva et al."1 d201 ll) and the emulator 



result, in Fig. [6] For z = 0, both ag ree at the 2% (low halo m ass) 
to 6% (high halo mass) level, the Bha ttacharva et al.l d201 ll) fit 
being slightly higher as expected. For z = 1, the discrepancy 
ranges from 4% to 17%. Here, overall the emulator estimate 
compared to the best-fi t power-law to the s i mulat ion result is 
slightly low, while the iBhattacharva et aL~1 d201 ll) fit slightly 
overestimates the simulation results. In other words, the actual 
simulation result lies in between the emulator pred iction and 
the fit. Overall, the agreement between the Bhattacharva et al. I 
d2011l) fit and the emulator is much better at z = 1 than the agree- 
ment between the emulator and the Bullock/Maccio fit. 

4.4. Cosmology Dependence of the c—M Relation 

Finally, we explore the sensitivity of the c-M relation to 
variations in cosmology. Since we now have a means of quickly 
and smoothly interpolating from one cosmology to another, we 
can simply vary the parameters that we incorporated into our 
design space in a regular grid. We divide each parameter range 
into five evenly spaced regions and vary only a single parame- 
ter at a time by keeping the other four parameters fixed at the 
midpoint of the parameter range. Our comparisons are always 
made with respect to this model at the midpoint, which has the 
parameters: n„,h 2 = 0.1375, fl b h 2 = 0.0225, n s = 0.95, w = -1 
and (78 = 0.75 8, and is subtracted from each c-M relation. Fig- 
ure [7] shows the results of this exercise at z = 0, with the val- 
ues of the cosmological parameters increasing as the shading 
increases from light to dark. The entire region is coloured to 
show the variation in the c-M relation across each parameter 
bin. The range of concentration variation changes as a func- 
tion of mass and the cosmological parameter being varied; the 
largest variation is of order unity. Unsurprisingly, the c-M re- 
lation increases with 0,„/i 2 and a% and decreases with w, which 
slows the rate of structure formation. We see little variation 
with fl b h 2 because the range of parameters allowed by the CMB 
constraints is already quite tight. Also, cluster-sized halos ap- 
pear to be relatively less sensitive to these changes in cosmol- 
ogy. Figure|7]also shows a clear degeneracy between ft m h 2 , erg, 
and n s in the c-M relation. 

5. CONCLUSION AND OUTLOOK 

In this paper, we have presented a new prediction scheme - in 
the form of an emulator - for the c-M relation for dark matter- 
dominated halos at the bright galaxy to cluster mass scales, cov- 
ering a range of 2 • 10 12 M© < M < 10 15 M© and a redshift range 
of z = to z = 1 . The emulator provides results for a large class 
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of wCDM cosmologies and is accurate at the ^5% level (bet- 
ter for lower redshifts, slightly worse for higher redshifts). The 
emulator enables consistent predictions to be made when test- 
ing for deviations from ACDM using clusters. This is partic- 
ularly important for cluster cosmology, since the behaviour of 
the c-M relation can vary by as much as 30% just by vary- 
ing the equation of state across the range — 1.3 < w < -0.7. By 
correctly including the cosmology dependence in the c-M rela- 
tion, the emulator improves on analytic modelling of halo pro- 
files, such as the 1-halo term used in the halo power spectrum. 
The performance of the emulator compares favourably with the 
other models for the c-M relation in the literature and outper- 
forms the Bullock/Maccio model across the redshift and mass 
range considered. 

Aside from predicting the mean c-M relation, the interesting 
and useful fact that across all 37 cosmologies considered, 1) the 
scatter in halo concentrations in individual mass bins is Gaus- 
sian, and 2) the corresponding standard deviation is given by 
roughly a third of the mean concentration value, means that the 
c-M emulator also includes within it the information regarding 
the concentration distribution at a given value of mass. 

The work in this paper is an example of how the cosmic 
calibration framework provides a means of estimating highly 
nonlinear quantities involving evolved structures from a lim- 
ited number of computationally expensive A/-body simulations. 
In the future, we will extend the number and range of cosmo- 
logical parameters to include more exotic phenomena, such as 
evolving dark energy, as a complement to upcoming dark en- 
ergy experiments. 
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